#ifndef Corr2_H
#define Corr2_H

#include <vector>
#include <string>

#include "Corr.h"
#include "CorrIO.h"
#include "dbg.h"
#include "ConfigFile.h"
#include "InputFile.h"
#include "BinData2.h"


template <int DC1, int DC2, int M>
class Corr2
{
public:
    Corr2(const std::vector<boost::shared_ptr<InputFile> >& files, const ConfigFile& params) : 
        _params(params), _vare(0.), _vare2(0.)
    {
        // Check for required parameters and set defaults for others.
        if (!params.keyExists("file_name")) myerror("Missing required parameter 'file_name'");
        if (DC1 != DC2 && !params.keyExists("file_name2")) 
            myerror("Missing required parameter 'file_name2'");

        std::string sep_unit_str = params.read("sep_units",std::string("arcsec"));
        double sep_units = (1. * ConvertUnitStr(sep_unit_str)) / angle::radians;

        if (!params.keyExists("nbins")) {
            // Then min_sep, max_sep, and bin_size are all required.
            if (!params.keyExists("min_sep")) myerror("Missing required parameter 'min_sep'");
            if (!params.keyExists("max_sep")) myerror("Missing required parameter 'max_sep'");
            if (!params.keyExists("bin_size")) myerror("Missing required parameter 'bin_size'");
            _minsep = double(params["min_sep"]) * sep_units;
            _maxsep = double(params["max_sep"]) * sep_units;
            _binsize = params["bin_size"];
            _nbins = int(ceil(log(_maxsep/_minsep)/_binsize));
        } else if (!params.keyExists("bin_size")) {
            // Then min_sep, max_sep, and nbins are all required.
            if (!params.keyExists("min_sep")) myerror("Missing required parameter 'min_sep'");
            if (!params.keyExists("max_sep")) myerror("Missing required parameter 'max_sep'");
            _minsep = double(params["min_sep"]) * sep_units;
            _maxsep = double(params["max_sep"]) * sep_units;
            _nbins = params["nbins"];
            _binsize = log(_maxsep/_minsep)/_nbins;
        } else if (!params.keyExists("max_sep")) {
            // Then min_sep, nbins, and bin_size are all required.
            if (!params.keyExists("min_sep")) myerror("Missing required parameter 'min_sep'");
            _minsep = double(params["min_sep"]) * sep_units;
            _nbins = params["nbins"];
            _binsize = params["bin_size"];
            _maxsep = exp(_nbins*_binsize) * _minsep;
        } else {
            // Then min_sep should not be specified.
            if (params.keyExists("min_sep")) 
                myerror("Only 3 of 'min_sep', 'max_sep', 'nbins', and 'bin_size should be "
                        "provided.");
            _maxsep = double(params["max_sep"]) * sep_units;
            _nbins = params["nbins"];
            _binsize = params["bin_size"];
            _minsep = _maxsep / exp(_nbins*_binsize);
        }
        dbg<<"nbins = "<<_nbins<<": min,maxsep = "<<_minsep<<','<<_maxsep<<std::endl;

        _logminsep = log(_minsep);
        _halfminsep = 0.5*_minsep;
        _minsepsq = _minsep*_minsep;
        _maxsepsq = _maxsep*_maxsep;

        double binslop = params.read("bin_slop",1.0);
        _b = binslop * _binsize;
        _bsq = _b * _b;

        _outputlevel = params.read("output_level",5);

        _data.resize(_nbins);

        // Create the field Cells from the input files.
        Assert(files.size() > 0);
        Assert(files.size() <= 2);
        std::vector<boost::shared_ptr<CellData<DC1,M> > > celldata;
        files[0]->buildCellData(_minsep,_binsize,celldata,_vare);
        _field1.reset(new Cell<DC1,M>(celldata));
        dbg<<"ngal = "<<getN()<<std::endl;
        dbg<<"vare = "<<_vare<<": sig_sn (per component) = "<<sqrt(_vare)<<std::endl;

        const int nfiles = files.size();
        if (nfiles > 2) myerror("nfiles > 2 not supported.");
        if (nfiles == 2) {
            std::vector<boost::shared_ptr<CellData<DC2,M> > > celldata;
            files[1]->buildCellData(_minsep,_binsize,celldata,_vare2);
            _field2.reset(new Cell<DC2,M>(celldata));
            dbg<<"file2: ngal = "<<getN2()<<std::endl;
            dbg<<"       vare = "<<_vare2<<": sig_sn (per component) = "<<sqrt(_vare2)<<std::endl;
        }
    }

    int getN() const { return _field1.get() ? _field1->getData()->n : 0; }
    int getN2() const { return _field2.get() ? _field2->getData()->n : 0; }
    double getVarE() const { return _vare; }
    double getVarE2() const { return _vare2; }

    void process()
    {
        Assert(_field1.get());
        if (_field2.get()) {
            process11(*_field1, *_field2, 1);
        } else {
            Assert(DC1 == DC2);
            process2(*_field1, 1);
        }
        finalize();
    }

    void write()
    {
        std::string unit_str = _params.read("sep_units",std::string("arcsec"));
        double units = (1. * ConvertUnitStr(unit_str)) / angle::radians;
        xdbg<<"minsep = "<<_minsep<<", units = "<<units<<std::endl;
        xdbg<<"minsep/units = "<<_minsep/units<<std::endl;
        if (_params.keyExists("e2_file_name")) {
            std::string filename = _params["e2_file_name"];
            std::ofstream os(filename.c_str());
            double smoothscale = _params.read("smooth_scale",0.);
            WriteEE(os,_minsep/units,_binsize,smoothscale,_data); 
        }
        if (_params.keyExists("m2_file_name")) {
            std::string filename = _params["m2_file_name"];
            std::ofstream os(filename.c_str());
            WriteM2(os,_minsep/units,_binsize,_data);
        }
    }

private:

    void process2(const Cell<DC1,M>& c1, int level);
    void process11(const Cell<DC1,M>& c1, const Cell<DC2,M>& c2, int level);
    void finalize();

    void directProcess11(const Cell<DC1,M>& c1, const Cell<DC2,M>& c2, const double dsq);

    const ConfigFile& _params;
    double _minsep;
    double _maxsep;
    int _nbins;
    double _binsize;
    double _binslop;
    int _outputlevel;
    double _logminsep;
    double _halfminsep;
    double _minsepsq;
    double _maxsepsq;
    double _b;
    double _bsq;

    double _vare;
    double _vare2;

    std::string _sepunits;

    boost::shared_ptr<Cell<DC1,M> > _field1;
    boost::shared_ptr<Cell<DC2,M> > _field2;

    std::vector<BinData2<DC1,DC2> > _data;
};

template <int DC1, int DC2, int M>
void Corr2<DC1,DC2,M>::directProcess11(
    const Cell<DC1,M>& c1, const Cell<DC2,M>& c2, const double dsq)
{
    XAssert(dsq > _minsepsq);
    XAssert(dsq < _maxsepsq);
    XAssert(c1.getSize()+c2.getSize() < sqrt(dsq)*_b + 0.0001);

    const double logr = log(dsq)/2.;
    XAssert(logr >= _logminsep);

    const int k = int((logr - _logminsep)/_binsize);
    xxdbg<<"logr = "<<logr<<", k = "<<k<<std::endl;
    XAssert(k >= 0); 
    XAssert(k<int(_data.size()));

    BinData2<DC1,DC2>& bindata = _data[k];

    bindata.directProcess11(c1,c2,dsq,logr);
}

template <int DC1, int DC2, int M>
void Corr2<DC1,DC2,M>::finalize()
{
    // TODO: Account corectly for vare2
    for(int i=0;i<_nbins;++i) _data[i].finalize(_vare);
    if (M == Sphere) {
        // Then all the measured distances were really secants.
        // Convert them to great-circle distances.
        for(int i=0;i<_nbins;++i) {
            double r = exp(_data[i].meanlogr);
            r = 2.*asin(0.5*r);
            _data[i].meanlogr = log(r);
        }
    }
    std::string unit_str = _params.read("sep_units",std::string("arcsec"));
    double units = (1. * ConvertUnitStr(unit_str)) / angle::radians;
    double log_units = log(units);
    for(int i=0;i<_nbins;++i) {
        if (_data[i].npair == 0.) _data[i].meanlogr = _logminsep+(i+0.5)*_binsize; 
        _data[i].meanlogr -= log_units;
    }
}

template <int DC1, int DC2, int M>
void Corr2<DC1,DC2,M>::process11(const Cell<DC1,M>& c1, const Cell<DC2,M>& c2, int level)
{
    const double dsq = DistSq(c1.getData()->pos,c2.getData()->pos);
    const double s1ps2 = c1.getAllSize()+c2.getAllSize();

    if (dsq < _minsepsq) if (sqrt(dsq)+s1ps2 < _minsep) return;
    if (dsq > _maxsepsq) if (sqrt(dsq)-s1ps2 > _maxsep) return;

    std::ostream* tempdbg = dbgout;
    if (dbgout) {
        if (level > _outputlevel) dbgout = 0;
        xdbg<< std::string(level,'.')
            << "Start P11: "<<c1.getData()->pos<<" -- "<<c2.getData()->pos
            << "   N = "<<c1.getData()->n<<","<<c2.getData()->n<<" d = "<<sqrt(dsq)<<std::endl;
    }

    // See if need to split:
    bool split1=false, split2=false;
    CalcSplitSq(split1,split2,c1,c2,dsq,_bsq);

    if (split1) {
        if (split2) {
            if (!c1.getLeft()) {
                std::cerr<<"minsep = "<<_minsep<<", maxsep = "<<_maxsep<<std::endl;
                std::cerr<<"minsepsq = "<<_minsepsq<<", maxsepsq = "<<_maxsepsq<<std::endl;
                std::cerr<<"c1.Size = "<<c1.getSize()<<", c2.Size = "<<c2.getSize()<<std::endl;
                std::cerr<<"c1.SizeSq = "<<c1.getSizeSq()<<
                    ", c2.SizeSq = "<<c2.getSizeSq()<<std::endl;
                std::cerr<<"c1.N = "<<c1.getData()->n<<", c2.N = "<<c2.getData()->n<<std::endl;
                std::cerr<<"c1.Pos = "<<c1.getData()->pos;
                std::cerr<<", c2.Pos = "<<c2.getData()->pos<<std::endl;
                std::cerr<<"dsq = "<<dsq<<", s1ps2 = "<<s1ps2<<std::endl;
            }
            Assert(c1.getLeft());
            Assert(c1.getRight());
            Assert(c2.getLeft());
            Assert(c2.getRight());
            process11(*c1.getLeft(),*c2.getLeft(),level+1);
            process11(*c1.getLeft(),*c2.getRight(),level+1);
            process11(*c1.getRight(),*c2.getLeft(),level+1);
            process11(*c1.getRight(),*c2.getRight(),level+1);
        } else {
            Assert(c1.getLeft());
            Assert(c1.getRight());
            process11(*c1.getLeft(),c2,level+1);
            process11(*c1.getRight(),c2,level+1);
        }
    } else {
        if (split2) {
            Assert(c2.getLeft());
            Assert(c2.getRight());
            process11(c1,*c2.getLeft(),level+1);
            process11(c1,*c2.getRight(),level+1);
        } else if (dsq > _minsepsq && dsq < _maxsepsq) {
            XAssert(NoSplit(c1,c2,sqrt(dsq),_b));
            directProcess11(c1,c2,dsq);
        }
    }


    if (tempdbg) {
        xdbg<<std::string(level,'.')<<"Done P11\n";
        dbgout = tempdbg;
    }
}

template <int DC1, int DC2, int M>
void Corr2<DC1,DC2,M>::process2(const Cell<DC1,M>& c12, int level)
{
    if (c12.getSize() < _minsep/2.) return;

    std::ostream* tempdbg = dbgout;
    if (dbgout) {
        if (level > _outputlevel) dbgout = 0;
        dbg<<std::string(level,'.')<<"Start P2: size = "<<c12.getSize()<<
            ", center = "<<c12.getData()->pos<<"   N = "<<c12.getData()->n<<std::endl;
    }

    xdbg<<std::string(level,'.')<<"P2 ("<<c12.getSize()<<") call Left P2\n";
    Assert(c12.getLeft());
    Assert(c12.getRight());
    process2(*c12.getLeft(),level+1);

    xdbg<<std::string(level,'.')<<"P2 ("<<c12.getSize()<<") call Right P2\n";
    process2(*c12.getRight(),level+1);

    xdbg<<std::string(level,'.')<<"P2 ("<<c12.getSize()<<") call P11\n";
    process11(*c12.getLeft(),*c12.getRight(),level+1);

    if (tempdbg) {
        dbg<<std::string(level,'.')<<"Done P2: size = "<<c12.getSize()<<
            ", center = "<<c12.getData()->pos<<std::endl;
        dbgout = tempdbg;
    }
}

#endif
